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SUMMARY 

Numerical investigation of transonic turbulent flows separated by 
streamline curvature and shock wave - boundary layer interaction is 
presented. The free stream Mach numbers considered are 0.4, 0.5, 0.6, 0.7, 
0.8, 0.825, 0.85, 0.875, 0.90, and 0.925. In the numerical method, the 
conservation of mass equation is replaced by a pressure correction equation 
for compressible flows and thus incremental pressure is solved for instead 
of density. The turbulence is described by a multiple-time-scale turbulence 
model supplemented with a near-wall turbulence model. The present numerical 
results show that there exists a reversed flow region at all free stream 
Mach numbers considered whereas various k-c turbulence models fail to 
predict such a reversed flow region at low free stream Mach numbers. The 
numerical results also show that the size of the reversed flow region grows 
extensively due to the shock wave - turbulent boundary layer interaction as 
the free stream Mach number is increased. These numerical results show that 
the turbulence model can resolve the turbulence field subjected to extra 
strains caused by the streamline curvature and the shock wave - turbulent 
boundary layer interaction and that the numerical method yields a 
significantly accurate solution for the complex compressible turbulent 
flow. 

*Work funded by Space Act Agreement C-99066-G. 
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coefficient for axial velocity correction equation 
coefficient for radial velocity correction equation 
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constant coefficient (=0,09) 

wall damping function for eddy viscosity equation 
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turbulent kinetic energy (k=kp + k t ) 
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axial and radial coordinates 
radial distance from the wall 
wall coordinate (=u r y/y) 

energy transfer rate of turbulent kinetic energy 

dissipation rate of turbulent kinetic energy 

dissipation rate of turbulent kinetic energy 

dissipation rate inside the near-wall layer 

von Karman constant (=0.41) 
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INTRODUCTION 


The transonic flow over an ax i symmetric curved hill [1] has received 
considerable attention in recent years as a bench mark test case to assess 
the capability of numerical methods as well as turbulence models to be used 
as design/analysis tools for fluid machinery. The transonic flow is 
schematically shown in Fig. 1. The boundary layer flow approaching the 
curved hill Is subjected to an extra mean flow strain rate generated by the 
streamline curvature. The development of the viscous force on the wall 
depends on the extra mean flow strain rate. As the fluid particle travels 
along the wall, the mean momentum Is dissipated by the strong viscous force 
and the flow eventually separates. As the free stream Mach number is 
further increased, a supersonic pocket is formed in the top region of the 
curved hill. As the strength of the shock wave Is increased with the 
increasing free stream Mach number, the reversed flow region grows 
extensively due to the shock wave - boundary layer Interaction. In 
numerical calculations of the transonic flow, correct prediction of the 
flow depends on the capability of a numerical method to resolve the 
compressible flow field which includes a supersonic flow region and a low 
Mach number reversed flow region and the capability of a turbulence model 
to properly resolve the turbulence field subjected to extra strain rates 
caused by the streamline curvature and the shock wave -boundary layer 
interaction. In this paper, calculations of the transonic flow at various 
free stream Mach numbers are made using a newly developed numerical method 
[2] and a multiple - time - scale turbulence model (hereafter abbreviated as 
the M-S turbulence model) [3,4]. A number of turbulence models, ranging 
from algebraic turbulence models to two-equation turbulence models 
incorporating a streamline curvature correction method, have been tested 
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and/or proposed in [5-7], Varying degrees of success have been reported in 
these references. The present numerical results are compared with these 
numerical results as well as the measured data. 

The Navier - Stokes equation solvers based on the pressure correction 
methods, also known as SIMPLE algorithms [8,9], are mostly used to solve 
Incompressible flows the domain of which can be discretized by an 
orthogonal mesh. Due to their strongly convergent nature, pressure 
correction methods have been used extensively to solve complex turbulent 
flows including chemically reacting turbulent flows. In the numerical 
method used herein, the pressure correction method has been extended to 
solve incompressible as well as compressible flows with arbitrary, complex 
geometries. The compressible flow equations are mostly solved by 
approximate factorization methods and flux splitting methods. The 
Beam-Warming method [10] and the MacCormack method [11] are representatives 
of the approximate factorization methods and the Steiger -Warming method 
[12] is a representative of the flux- splitting methods. These methods were 
originally developed to solve the Euler equations and were then extended to 
include the viscous term to solve the Navier- Stokes equations. A few 
differences exist between the two classes of methods. In the latter class 
of methods, the density is solved for as a primary variable and the 
pressure is obtained from the equation state. For incompressible flows, the 
pressure no longer depends on the density and hence the latter class of 
methods fails for incompressible flows. These methods can also be extended 
to solve Incompressible flows by including an artificial compressibility 
into the governing flow equations [13]. On the other hand, in the pressure 
correction methods, the Incremental pressure is solved for as a primary 
variable, hence the method is valid for both incompressible and 


5 



compressible flows. Another difference between the two classes of methods 
can be found in the way the second order diffusion term is treated. In the 
pressure correction methods, the diffusion term is incorporated into the 
stiffness matrix while, in the other class of methods, the diffusion term 
is incorporated into the system of equations as the load vector term. For 
turbulent flows with extensive recirculation zones, the pressure correction 
methods may be numerically more stable than the other class of methods, 
conceptually; however, the pressure correction methods have mostly been 
used for incompressible flows and the approximate factorization methods and 
the flux splitting methods have mostly been used for compressible flows. 
Therefore, definitive advantages and disadvantages of these two classes of 
methods can not be discussed with confidence as yet. 

A few papers to extend the SIMPLE method to solve compressible flows 
with complex geometries have appeared in recent years [14-16]. Some 
difficulties have been encountered in the course of these studies. One 
difficulty was identifying a suitable grid layout to solve the 
Navier- Stokes equations defined on complex geometries. In [14], a 
collocated grid layout was used and an artificial dissipation was included 
to prevent velocity-pressure decoupling. In the present numerical method 
and in [16], the velocities are located at the same grid points and the 
pressure is located at the centroid of pressure control volume formed by 
the four adjacent velocity grid points. This grid layout may yield a 
velocity-pressure decoupled solution if used together with the standard 
pressure correction procedure [17]. The mechanism that leads to the 
velocity-pressure decoupled solution is heuristically shown in [18], In 
[14], the velocity-pressure decoupling was eliminated by using a 
non-conforming control-volume for mass imbalance calculation. In this case, 
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an uncertainty caused by the use of a non- conforming domain for the mass 
imbalance calculation needs to be further investigated. In the present 
numerical method, the velocity-pressure decoupling is eliminated by 
treating the pressure correction equation as a standard partial 
differential equation rather than treating it as a constraint condition. 
Further details are discussed in the following section. Another difficulty 
was to find a simple, yet strongly convergent, pressure correction equation 
valid for both incompressible and compressible flows. The capability to 
solve compressible flows is achieved by Including a convective pressure 
correction term into the disturbed conservation of mass equation in one 
form or another. However, the numerical procedures to solve each 
compressible form pressure correction equation differ significantly from 
each other. A multi-step pressure correction algorithm was used in [14] , 
and the SIMPLE-R [8] and the SIMPLE-C [19] were used in [15] and [16], 
respectively. In these methods, the density was also corrected from the 
incremental pressure. In the present method, only the pressure and the 
velocity are corrected from the incremental pressure and the density is 
obtained from the equation of state. Thus the present method Is simpler 
than the other methods, and the multi-step pressure correction algorithm 
[14] is more involved than the other methods considered herein. The 
accuracy and the convergence nature of the present numerical method has 
been tested by solving a number of example flows. The example flows 
considered In [2] include: a developing channel flow, a developing pipe 
flow, a two-dimensional laminar flow in a 90 degree bent channel, polar 
cavity flows, and a turbulent supersonic flow over a compression ramp. More 
calculations of various complex turbulent flows using the same numerical 
method can be found in [20-22]. It can be seen from these numerical results 
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that the present numerical method yields accurate computational results 
even when highly skewed, unequally spaced, curved grids are used. 

It has long been known that the turbulent transport is related to the 
time scale of energy containing large eddies and the dissipation of 
turbulent kinetic energy is related to the time scale of fine scale eddies 
in the dissipation range [23]. In M-S turbulence models, the turbulent 
transport of mass and momentum Is described using the time scale of the 
large eddies and the dissipation rate is described using the time scale of 
the fine-scale eddies. Due to the physically consistent nature of the M-S 
turbulence models, these turbulence models are expected to yield 
significantly improved computational results compared with the 
single- time- scale turbulence models. However, the first M-S turbulence 
model [24] did not quite come up to the expectations due to a few 
shortcomings in the closure model. These shortcomings and a few differences 
between the two M-S turbulence models are discussed in the following 
section for the record. On the other hand, the present M-S turbulence model 
yields significantly improved computational results than the 
single- time - scale turbulence models for a number of complex turbulent flows 
[3,20,21]. These complex turbulent flows include: a wall-jet, a 
wake-boundary layer interaction flow, a turbulent flow over a 
backward- facing step, a confined coaxial swirling jet, turbulent flows over 
a strongly curved surface, and reattaching shear layers In a divergent 
channel. Calculation of more complex turbulent flows are in progress. 

In numerical calculations of turbulent flows, the near-wall turbulence 
is usually described using the wall functions [25], two- or multi-layer 
turbulence models [26,27], and low Reynolds number turbulence models [28], 
In the present study, the near-wall turbulence is described by a “partially 
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low Reynolds number approach" [4]. In this near-wall turbulence model, only 
the turbulent kinetic energy equations are extended to include the 
near -wall low turbulence region and the energy transfer rate and the 
dissipation rate inside the near-wall layer are obtained from algebraic 
equations. The algebraic equations were obtained from a k-equation 
turbulence model [29]. The advantages of the present near-wall turbulence 
model over the low Reynolds number turbulence models can be described as 
follows. The turbulence length scale of boundary layer flows is strongly 
related to the normal distance from the wall. This characteristic of the 
wall bounded turbulent flows can be described quite naturally by empirical 
algebraic equations. The low Reynolds number turbulence models can also be 
used to describe the wall bounded turbulent flows; however, more grid 
points have to be used to resolve the steep dissipation rate in the 
near-wall region. More detailed discussion on the advantages and 
disadvantages of various near-wall turbulence models, the development of 
the present near-wall turbulence model, and its application to fully 
developed turbulent channel and pipe flows can be found in [4], It is also 
shown in [4] that the near-wall turbulence model can resolve the over- shoot 
phenomena of the turbulent kinetic energy and the dissipation rate in the 
region very close to the wall. Incorporation of the near-wall turbulence 
model into a k-£ turbulence model and its application to complex turbulent 
flows such as a supersonic turbulent flow over a compression ramp and a 
transonic flow over an axlsymmetric curved hill can be found in [2] and 
[22] , respectively. 

REYNOLDS AVERAGED NAVIER- STOKES EQUATIONS AND NUMERICAL METHOD 
The compressible turbulent flow equations are given as; 
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and the density is obtained from the perfect gas law given as p=/)RT. A 
turbulent Prandtl number (cr^) of 0.75 was used for the energy equation. The 
molecular viscosity and the thermal conductivity were obtained from 
Sutherland's laws given as [30]; 
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where « 1.716 x 10 Kg/m-sec, T Q — 273.1° Kelvin, S — 110.6° Kelvin; 
and 
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where k Q - 0.0264 Kg/m-K, T 0 - 273.1° Kelvin, and S - 194.4° Kelvin. 

The specific heat was obtained from a curve-fitted 4-th order polynomial, 
see [31] for details. 

The pressure correction equation valid for both incompressible and 
compressible flows is described below. As in the standard pressure 
correction method, the density, the velocity, and the pressure are 
decomposed as; 
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Substituting eqs . (7-9) into eq. (1) yields: 

V- (p'V*) + V- (pV ) + V(p'V') V- (p*V*) (11) 

The third term on the left hand side of eq. (11) is neglected for 
simplicity in any of the pressure correction algorithms discussed below. 
The incremental pressure is related to the incremental density and the 
incremental velocities as; 
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where eq. (12) is obtained from the equation of state and eqs. (13-14) are 
obtained from the discrete u- and v-momentum equations, respectively. 
Substituting eqs. (12-14) into (11) yields, after some rearrangement; 


a 

r * ■> 

u 

i a 

r * > 

V 

a 

a P 'l 

i a 

+ 3p r 


— P' 

+ 

r— p' 

— 

P A u 

— 

rpA v 

dx 

l RT J 

r 3r 

t RT 

dx 

3x J 

r 3r 

dr J 


12 


In the present numerical method, all flow variables, except pressure, 
are located at the same grid points and the pressure is located at the 
centroid of the cell formed by the four neighboring grid points. A few 
remarks on the pressure correction algorithm are in order for clarity. In 
the more standard pressure correction algorithms, the discrete pressure 
correction equation Is obtained by directly substituting the discrete form 
incremental pressure - incremental velocity relations, eqs . (13-14), into 
eq. (11). In this case, the discrete pressure correction equation for a 
pressure grid point is given as a nine-diagonal system of equations for 
rectangular grids. This pressure correction equation can yield a 
velocity-pressure decoupled solution as discussed in [17,18]. Also this 
pressure correction equation is not diagonally dominant. On the other hand, 
in the present method, the continuous form pressure correction equation, 
given as eq. (15), is solved for incremental pressure. In this case, the 
discrete pressure correction equation is given as a five - diagonal system of 
equations for rectangular grids. This discrete pressure correction equation 
is strongly diagonally dominant even for highly skewed grids. Even the 
slightest symptom of velocity-pressure decoupling Is not observed with the 
present pressure correction algorithm. 

The capability to solve compressible flows with shock waves is 
achieved by the convective incremental pressure terms, the first two terms 
in the left hand side of eq . (15). These two terms properly take Into 
account the hyperbolic nature of supersonic flows, and enable the capture 
of shock waves. In low Mach number flows and in the near-wall boundary 
layers of supersonic flows, the variation of density mostly depends on the 
local temperature. However, the dependence of density on temperature 
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has been ignored in deriving the convective terms for simplicity, see eq . 
(12). Yet, the dependence of density on temperature is clearly resolved in 
the present numerical method, since the incremental pressure is driven only 
by the mass imbalance evaluated from the conservation of mass, the right 
hand side of eq . (15), and the density is obtained from the equation of 
state. In fact, it can be seen in [32] that low Mach number flows can be 
solved even without the convective incremental pressure terms. Equally 
importantly, use of the simplified incremental pressure -incremental 
density relationship still yields a rapidly convergent solution as shown in 
the following section. Also note that the present pressure correction 
algorithm is significantly simpler than the multi-step pressure correction 
algorithm [14]. 

In solving the system of equations, the power- law upwinding [8] is 
used for all convection-diffusion equations except for the pressure 
correction equation. The upwind differencing [8] is used for the pressure 
correction equation. In the region very close to the wall, highly fine 
grids need to be used to resolve the thin boundary layer properly. In this 
region, the numerical method is second order accurate; however, the method 
becomes first order accurate in the free stream region where the mesh is 
coarse. Each differential equation is solved sequentially until the 
relative error for each flow variable becomes smaller than the prescribed 
convergence criterion and the mass imbalance in eq. (15) becomes 
negligible. The incremental pressure is obtained by solving eq. (15). In 
solving the discrete system of equations, the off-diagonal terms may be 
moved to the load vector term and the resulting system of equations can be 
solved using a tri-diagonal matrix algorithm (TDMA) . The corresponding 
incremental velocities are obtained from eqs . (13-14). The flow variables 
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are updated using eqs . (8-10), and these updated flow variables are used in 

computing the new current flow variables by solving eqs. (2-4) and the 
turbulence equations. 


TURBULENCE EQUATIONS 

The M-S turbulence model supplemented with the near-wall turbulence 
model is summarized below for completeness. The turbulent kinetic energy 
and the energy transfer rate equations for the energy containing large 
eddies are given as; 
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where Pr^i^S* is the production rate. The turbulent kinetic energy equation 
and the dissipation rate equation for the fine scale eddies are given as: 
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The eddy viscosity is given as; 


k 2 

M t “ P c fif — 
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( 20 ) 


The turbulent kinetic energy equations, eqs . (16) and (18), are defined for 
the entire flow domain while the energy transfer rate equation, the 
dissipation rate equation and the eddy viscosity equation are valid for the 
flow domain away from the near-wall region. The turbulence model constants 
are given as; <7^-0. 75, a kt =( ^-^» Cp]_=0.21, Cp2~1.24, 

c p3 “ 1.84, c tl =0 . 29 , c t2 “ 1.28, and c t3 -1.66. Details on the present M-S 

turbulence model can be found in [3], 

The energy transfer rate and the dissipation rate inside the near-wall 

layer are given as; 
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and the eddy viscosity for the near-wall layer is given as; 


k 2 

"t “ c /if f M — 
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( 23 ) 


where f^-1-1 . /exp (Aj^R,- + A 2 R t 2 ). The coefficients A x and A 2 are given as 
0.025 and 0.00001, respectively, see [4]. 

The domain for each differential equation is shown schematically in 
Figure 2. For wall bounded turbulent flows, the equilibrium region extends 
from y + ~30 to y + «300. Thus the partition between the near-wall region and 
the fully turbulent outer region can be located between y + greater than 30 
and less than 300 approximately. Recall that the present near-wall 
turbulence model and the k- equation turbulence model [29] from which the 
present near-wall turbulence model is derived are valid for the entire flow 
domain of equilibrium boundary layer flows. Thus the computational results 
do not depend appreciably on the location of the partition. However, if the 
partition is located too far away from the wall (i.e., y + >1000) , then the 
numerical results in the near-wall region may become similar to those 
obtained using a k-equation turbulence model. 

A few differences between the present M-S turbulence model and that of 
[24] are summarized in this paper for the record. Firstly, the eddy 
viscosity equation in [24] is given as; 



(24) 


Numerical calculations of complex turbulent flows showed that the ratio of 
k t /kp can vary significantly in regions where the turbulence is in a 
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strongly Inequi librium state. Eq. (24) implies that the small scale eddies 
contained in the dissipation range may not contribute significantly to the 
turbulent transport of mass and momentum. This anomaly can be cured if the 
partition between the large eddies and the small eddies is moved toward the 
very high wave number region so that the k t may be negligibly small in all 
occasions. However, in this case, the multiple- time -scale turbulence model 
can be reduced to a single- time-scale turbulence model as discussed in 
[24]. Secondly, in the present M-S turbulence model, the variable energy 
transfer functions were obtained from a physical dimensional analysis [30]. 
On the other hand, the M-S turbulence model in [24] contains such a 
variable energy transfer function only in the energy transfer rate 
equation. Hence the load functions of the energy transfer rate equation and 
the dissipation rate equation lack symmetry. Thirdly, in the present M-S 
turbulence model , the model constants were obtained by solving a five by 
five system of equations obtained by transforming the M-S turbulence 

equations into asymptotic equations for the decay rate of grid turbulence 

[33] and the growth rate of turbulence intensity [34], Lastly, of practical 
importance, the eddy viscosity equation given as eq. (24) is inconsistent 
with the near wall analysis unless k^ vanishes in the near-wall equilibrium 
region, see [3]. In application to complex turbulent flows, arbitrary 
ratios of k t /kp were used as a near wall boundary condition together with 
the standard wall functions [24,35]. A wall function for the M-S turbulence 
model obtained from a near-wall analysis is given in [3], if any wall 

function need to be used. Also an arbitrary ratio of k t /k p was used as an 

inlet boundary condition in a number of boundary layer calculations [35]. 

In this case, the calculated shear layer expands rapidly so that the 
turbulence field can adjust itself to the ill -posed inlet boundary 
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condition, see [3]. 


COMPUTATIONAL RESULTS 

The measured data for the transonic flow over an axisymmetr ic curved 
hill at various free stream Mach numbers can be found in [1,5,36]. In the 
experiment, an axisymmetric circular-arc bump of thickness 1.9 cm and a 
chord length of 20.3 cm was attached 60 cm downstream of a circular 
cylinder with the external diameter of 15.2 cm. The numerical results for 
M^-0.875 and R e -13 . 2xl0 6 /m are compared with the measured data presented in 

[36] . The boundary layer thickness of the approaching transonic flow was 
0.01 meters for M^-0.875. The other numerical results for the rest of the 
free stream Mach numbers at R e —10xl0^/m are compared with the measured data 
given in [ 5 ] . 

In the following calculations , the inlet boundary is located at one 
chord length upstream of the forward corner of the bump; and the exit 
boundary, at one chord length downstream of the rear end of the bump. Some 
degree of uncertainty that may be caused by numerical diffusion and 
inadequate grid size can always exist in any numerical analysis. To reduce 
the numerical uncertainty, three different meshes (78x53, 108x65, and 
145x65 grid points in the axial and radial directions, respectively) have 
been used in the present study. The computational results obtained using 
the first two grids differed by no more than a few percent; and the latter 
two, by no more than one percent. The computational results presented 
herein were obtained using the finest grid shown in Figure 3. The Inlet 
boundary condition for the axial velocity and the turbulent kinetic energy 
were obtained from experimental data for a fully developed flat plate flow 

[37] . The non-dimensional velocity and the turbulent kinetic energy 
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profiles were scaled to yield a boundary layer thickness of 0,01 meters at 
the inlet boundary. Uniform static pressure and uniform enthalpy were also 
prescribed at the inlet boundary. The no- slip boundary condition for 
velocities, a vanishing turbulent kinetic energy, and a constant 
temperature which corresponds to the free stream stagnation temperature 
were prescribed at the solid wall boundary. The free stream flow condition 
was prescribed at the top boundary, and a vanishing gradient boundary 
condition was used for all flow variables at the exit boundary. The 
partition between the near-wall layer and the external region was located 
at approximately one percent of the boundary layer thickness away from the 
wall. Thus the partition is located at y + ~80 (for M^O.875) at the inlet 
boundary and 11 grid points were allocated inside the near-wall layer. The 
mesh size of the first grid point on the bottom wall was Ay + ~1,25 and the 
grid size in the normal to the wall direction was increased by a factor of 
approximately 1.15. The initial guess was obtained by extending the inlet 
boundary condition in the axial direction. An uncertainty that can be 
caused by the location of the inlet boundary and the inlet boundary 
conditions is clarified by comparing the numerical results with the 
measured data at x/c--0.25. 

The convergence history for M^-0.875 is shown in Figure 4. Each 
iteration consists of 7 sweeps of the pressure correction equation and 3 
sweeps for the rest of the flow equations in the axial and in the radial 
directions, respectively. The pressure was updated using an 
under- relaxation factor of 0.57; and the rest of the flow variables, using 
an under- relaxation factor of 0.47. The relative error for each flow 
variable shown in Figure 4(a) is defined as; 
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Ri - Maxl|(a n+1 - a n )/A n+1 |. J-l.N) 

* j ^ » j ^ 


( 25 ) 


where the superscript n denotes the iteration level; the subscript 
i={u,v,p,T) denotes each flow variable; the subscript j denotes each grid 
point; N denotes the total number of degrees of freedom for each flow 
variable; and denotes the maximum magnitude of the i-th flow variable. 
The convergence histories for the other scalar variables (kp, e pi k t , and 
e t ) are almost the same as that of temperature. The mass imbalance shown in 
Figure 4(b) is defined as; 

N c 

R 2 ” VN C [ l (V-( P V)) 2 C ] (26) 

c-1 

where N c is the total number of the pressure control volumes. The 
"practically” converged solution was obtained in approximately 1000 
iterations for all the free stream Mach numbers considered. Note that, in 
control -volume based finite difference methods, the discrete system of 
equations is derived by integrating the governing differential equations 
over each control volume. For curvilinear grids, the required number of 
interpolations to obtain flow variables at the cell boundaries is 
significantly reduced by using the present grid layout than the one used in 
[18]. The strongly convergent nature of the present numerical method is 
partly attributed to the grid layout which required fewer interpolations 
and the pressure correction algorithm which yields a diagonally dominant 
system of equations even when highly skewed meshes are used. It can also be 
found in Figure 4(b) that the mass imbalance converges almost 
monotonically . Such a monotonically convergent nature is attributed to the 
use of the under- relaxation . With the use of these under-relaxation 
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parameters, divergence or convergence to an erroneous solution was not 
encountered. 

The calculated iso-Mach lines are shown in Figure 5, where the 
incremental Mach number between the contour lines is constant for each free 
stream Mach number. It can be seen in this figure that a small supersonic 
pocket first appears at M^-0.80. The iso-Mach lines are almost symmetric at 
low free stream Mach numbers and the symmetry is lost as the size of the 
supersonic pocket grows with the increasing free stream Mach number. It can 
be seen from these figures that the present numerical method can cleanly 
resolve the transonic flows from the low to the high transonic free stream 
Mach number. The size of the supersonic pocket for M^-0.925 also compares 
favorably with that obtained using the MacCormack scheme [5] . The 
calculated shock is slightly more spread out than that of [5] since the 
present numerical method becomes first order accurate in the free stream 
region where coarse grids are used. However, the slightly smeared shock 
does not impair the numerical results appreciably as can be seen in this 
section. 

The calculated static pressure distributions on the wall are compared 
with the measured data as well as the other numerical results in Figure 6. 
It can be seen In the figure that the pressure distributions on the wall 
obtained by the MacCormack scheme using the King-Johnson algebraic 
turbulence model [5] (hereafter abbreviated as K-J model) compare most 
favorably with the measured data. The present numerical results show that 
the calculated shocks are located somewhat downstream of the measured data 
for all free stream Mach numbers. This discrepancy is attributed to the 
turbulence model which slightly under-estimates the Reynolds stress. For 
M^-0.875, the present result compares more favorably with the measured data 
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than does the one obtained by the MacCormack scheme using the 
Wilcox-Rubesin Turbulence model (hereafter, abbreviated as W-R model) [7], 
However, it can be found in [6] that the W-R model yields substantially 
improved numerical results for the wall pressure and the mean velocity if 
the free stream condition is prescribed at the outer boundary* 

The streamline and the pressure contour lines for M^-0.875 are shown 
in Figures 7(a) and 7(b), respectively. These contour lines also indicate 
that the present numerical method can cleanly resolve the transonic 
turbulent flow. 

The calculated separation and reattachment locations are compared with 
the measured data and the other numerical results in Figure 8. It is shown 
in this figure that the present method successfully predicts the existence 
of the reversed flow region at all free stream Mach numbers. At low free 
stream Mach numbers, the present results compare more favorably with the 
measured data than do those obtained by the MacCormack scheme using the K-J 
turbulence model [5]. As the free stream Mach number is increased, the 
present method slightly under-predicts the size of the reversed flow region 
compared with the measured data and the numerical results obtained using 
the K-J turbulence model. This under-prediction of the reversed flow region 
is a result of the calculated shocks located slightly downstream of the 
measured data. It is also shown in this figure that the Jones -Launder k-e 
turbulence model [5] and a k-e turbulence model supplemented with a 
streamline curvature correction method [6] fail to predict the reversed 
flow region at low free stream Mach numbers. These turbulence models also 
under-predict the reattachment locations at high free stream Mach numbers 
as shown in the figure. 

The mean velocity profiles for M^-0.875 at five axial locations are 
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compared with the measured data as well as the other numerical results in 
Figure 9. At x/c--0.25, the calculated mean velocity profile and the 
measured data compare favorably with each other, which indicates that the 
inlet boundary condition used in the present study is a good approximation 
to the experiment at the inlet boundary. It can be seen in the figure that 
all numerical results exhibit fair comparison with the measured data. In 
particular, the mean velocity profile obtained using the W-R turbulence 
model [7] compares less favorably with the measured data than the other 
numerical results. Again, it can be found In [6] that the W-R model yields 
an improved mean velocity profile if the free stream condition is 
prescribed at the outer boundary. It is also shown in the figure that a k-e 
turbulence model incorporating an improved wall function [25] 
under-predicts the magnitude of mean velocity at x/c=Q,75 and 0.875 [39]. 
This under-prediction in the mean velocity is caused by the over-predicted 
Reynolds stress at the same axial locations. 

The Reynolds stress profiles for M^-0.875 at five axial locations are 
shown in Figure 10, It can be seen in the figure that the calculated and 
the measured Reynolds stress profiles at x/c-*-0.25 compare favorably with 
each other, which, again, indicates that the inlet boundary condition used 
in the present study Is a good approximation to the experiment at the inlet 
boundary. At low free stream Mach numbers for which the shock wave 
-boundary layer interaction do not exist, the flow separation is caused by 
the turbulent shear stress developing over the forward part of the curved 
hill [38]. A successful prediction of such a flow depends on the capability 
of a turbulence model to correctly describe the turbulence field subjected 
to the streamline curvature [20], As shown in this figure, the present 
numerical results compare more favorably with the measured Reynolds stress 
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than the other numerical results at x/c-0.69 and 0.75. However, the 
magnitude of the present numerical results is slightly smaller than the 
measured data at these locations. The calculated shocks and separation 
points are located slightly downstream of the measured data at high free 
stream Mach numbers, which is attributed to the slightly under-predicted 
Reynolds stress in the forward part of the curved hill. It Is also shown in 
the figure that the k-e turbulence model with an Improved wall function 
[39] significantly over-estimates the Reynolds stress at x/c=0.75. Inside 
the reversed flow region, x/c— 1.0, the present numerical result compares 
less favorably with the measured data than does the one obtained using the 
K-J turbulence model. This under-prediction in the magnitude of the 
Reynolds stress Is attributed to the calculated shock and the separation 
point which are located slightly downstream of the measured data. It is 
also shown in the figure that the Reynolds stress profile at x/c=1.38 
obtained using a modified K-J turbulence model, with the modifications 
restricted to the reversed flow region, compares more favorably with the 
measured data than other numerical results. However, this model is shown to 
over-predict the Reynolds stress near the outer edge of the reversed flow 
region. 

In the course of the development of the present numerical method, the 
same transonic flow at M CO =0.875 has been solved using a k-e turbulence 
model supplemented with the same near-wall turbulence model as used in the 
present study. The k-e turbulence model [22] also successfully predicted 
the existence of the reversed flow region at all free stream Mach numbers. 
However, the size of the reversed flow region for each free stream Mach 
number was approximately 10 percent smaller than the present result. The 
smaller size of the reversed flow region is due to the incapability of k-c 
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turbulence models to resolve turbulence fields in strongly inequilibrium 
state such as that inside the reversed flow region. The turbulent kinetic 
energy obtained using the present M-S turbulence model is almost the same 
as that obtained using the k-e turbulence model. Both turbulence models 
significantly under-predicted the turbulent kinetic energy, see [22] for 
more details. 


CONCLUSIONS 

Calculations of turbulent transonic flows with a control -volume method 
based on a pressure correction method were presented. The turbulence was 
described by a multiple- time -scale turbulence model supplemented with a 
"partially low Reynolds number" near-wall turbulence model. 

The numerical results showed that the supersonic pocket on the top of 
the curved hill first appeared at the free stream Mach number of 0.80 and 
that the supersonic pocket became larger as the free stream Mach number was 
further increased. The numerical results also showed that there exists a 
reversed flow region at low free stream Mach numbers and that the size of 
the reversed flow region grows extensively as the free stream Mach number 
is increased. Thus the numerical method was shown to yield a significantly 
accurate solution for the complex compressible turbulent flow including the 
supersonic pocket and the nearly Incompressible low Mach number reversed 
flow region. 

For turbulent flows over a curved hill, the mean flow is subjected to 
extra strains caused by the streamline curvature. The development of the 
turbulence field over such a curved surface mostly depends on the extra 
strains. The capability to predict the reversed flow region in turbulent 
flows over a curved hill rests on the capability of a turbulence model to 
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properly resolve the turbulence field. It was shown that the present 
turbulence model can predict the reversed flow region at low free stream 
Mach numbers while the Jones -Launder k-c turbulence model and a k-€ 
turbulence model supplemented with a streamline curvature correction method 
fail to predict the reversed flow region [5,6]. The present numerical 
results also show the extensive growth of the reversed flow region caused 
by the shock wave - turbulent boundary layer interaction at high free 
stream Mach numbers. These numerical result compare favorably with the 
measured data and the other numerical results obtained using the 
King- Johnson turbulence model [5]. 
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FIGURE 2. - DOMAIN FOR EACH PARTIAL DIFFERENTIAL EQUATION. 



FIGURE 3. - DISCRETIZATION OF THE FLOW DOMAIN <1^ X 65 P€SH). 
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FIGURE 8. - FLOW SEPARATION AND REATTACHMENT iOCATIONS. 
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FIGURE 10. - REYNOLDS STTESS PROFILES. 
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